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Jonas Vlictinck,^ Jan Ryckcbusch/ and Kris Van Houckc^ 
' Department of Physics and Astronomy, Ghent University, Proeftuinstraat 86, 9000 Gent, Belgium 

We report on a study of a spin-down impurity strongly coupled to a spin-up Fermi sea (a so- 
called Fermi polaron) with the diagrammatic Monte-Carlo (DiagMC) technique. Conditions of zero 
temperature and three dimensions are considered for an ultracold atomic gas with resonant inter- 
actions in the zero-range limit. A Feynman diagrammatic series is developed for the one-body and 
two-body propagators providing information about the polaron and molecule channel respectively. 
The DiagMC technique allows us to reach diagram orders that are high enough for extrapolation to 
infinite order. The robustness of the extracted results is examined by checking various resumma- 
tion techniques and by running the simulations with various choices for the propagators and vertex 
functions. It turns out that dressing the lines in the diagrams as much as possible is not always 
the optimal choice. We also identify classes of dominant diagrams for the one-body and two-body 
self-energy in the region of strong interaction. These dominant diagrams turn out to be the leading 
processes of the strong-coupling limit. The quasiparticle energies and Z-factor are obtained as a 
function of the interaction strength. We find that the DiagMC results for the molecule and polaron 
properties are very similar to those obtained with a variational ansatz. Surprisingly, this variational 
ansatz gives very good predictions for the quasiparticle residue even when this residue is significantly 
smaller than one. 
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I. INTRODUCTION 

The notion of a 'bare' partiele loses its significance once 
it is strongly coupled to a medium. Landau introduced 
the notion of a quasiparticle whose properties may be 
very different from those of a bare particle^. The most 
prominent example is an electron moving in a crystal: the 
electron displaces the nearby ions and carries this distor- 
tion with it. The presence of the phonon cloud changes 
the mass and energy of the electron, that is dubbed as 
'polaron—. More generally, a polaron arises whenever a 
quantum impurity is strongly coupled to an environment. 
These quantum-mechanical quasiparticlcs play a key role 
in the low-energy behavior of a macroscopic quantum liq- 
uid. 

In recent years, the field of ultracold atoms has pro- 
vided an exciting framework for studying polaronic ef- 
fects. A key idea is that models designed for describing 
the rich and non-trivial structure of the solid state, can 
be emulated in a clean and controllable manner with ul- 
tracold atoms. For example, so-called Fermi polarons^"— , 
spin-down impurities that are strongly coupled to a 
spin-up Fermi sea, can be created in a degenerate two- 
component atomic Fermi gas when going to the limit of 
strong spin-imbalance close to a Feshbach resonance. The 
impurity is coherently dressed with particle-hole excita- 
tions of the Fermi sea. The properties of the Fermi po- 
laron are important for the quantitative understanding 
of a strongly imbalanced Fermi gas^. 

In this paper, we focus on the 'attractive Fermi po- 
laron', with an attractive interaction between the impu- 
rity and the fermions of the bath. A recent experiment 
using an ultracold gas of ^Li atoms in three dimensions 
revealed the existence of Fermi polarons through a nar- 
row quasiparticle peak in the impurities' radio-frequency 



(rf) spectrum^. At a critical interaction strength, the 
disappearance of this peak was interpreted as a transi- 
tion from polaronic to molecular binding, when the im- 
purity and an atom of the sea form a two-body bound 
state. Such a transition had theoretically been predicted 
in three dimensions (3D) by Prokof'ev and Svistuno^*^. 
To determine the transition point, they developed a dia- 
grammatic Monte Carlo technique (DiagMC) capable of 
solving the Fermi polaron modcl^-i^. Calculations of the 
ground-state energy showed that for a sufficiently strong 
attraction between the impurity atom and the atoms of 
the spin-up Fermi sea, a molecular state becomes ener- 
getically favorable. The crossing point was found at an 
interaction strength {kpa)c = 1.11(2), with kp the Fermi 
momentum of the spin-up sea and a the s-wave scatter- 
ing length. A variational treatment developed by Chevy 
based on an expansion up to single particle-hole excita- 
tions on top of the unperturbed Fermi sea turned out 
to be remarkably accurate^. A combination of Chevy's 
ansatz with a variational wave function in the molecular 
limitiSi"— also revealed the polaron-to-molecule transi- 
tion, very close to the DiagMC result. 

In the present work, we study the quasiparticle prop- 
erties of the Fermi polaron problem in 3D with the Di- 
agMC technique^i^. This technique evaluates a series of 
Feynman diagrams for the one-particle and two-particle 
proper self- energies. A full description of the DiagMC 
algorithm was presented in Ref4. Building on the work 
of Rcf4 we have implemented the DiagMC algorithm in- 
dependently. We explore various DiagMC schemes^ and 
series resummation methods to check the robustness of 
the results against the possible uncertainties of summing 
the series. First, we confirm the transition point. Next, 
we calculate the quasiparticle residue which we compare 
to experimental data and variational results. The quasi- 



2 



particle residue, or Z-factor, gives the overlap of the non- 
interacting wave function and the fully interacting one, 

Z, = |(<^|0,,F5(7V^))P, (1) 

with l*?^^) the fully interacting ground state and 
IO4,, FS{N-f)) a free spin-down atom carrying momentum 
p = in a non-interacting Fermi sea FS of TV-j- spin- 
up atoms. The spin-up atoms are non-interacting since 
p-wave scattering is negligible. The residue reflects the 
impurity's probability of free propagation. 

The outline of the paper is as follows. In Section llll we 
introduce the model and the structure of the Feynman 
diagrammatic expansion. In Section IIIII we discuss the 
results of the numerical calculations. Thereby, we inves- 
tigate on how the results depend on the choices made 
with regard to the diagrammatic series, like the use of 
bare versus dressed propagators. Also the resummation 
of the diagrammatic series is discussed in depth. The 
results for the quasiparticle properties, like the residue, 
are the subject of Section IIVI 



II. MODEL AND DIAGRAMMATIC 
STRUCTURE 

We consider a dilute two-component gas of ultra- 
cold fermionic atoms interacting via the van der Waals- 
potential. The Hamiltonian has a kinetic and interaction 
term 

k,o-=ti 
k,k',q 

(2) 

The operators cj^.^. (C]^.^) create (annihilate) fermions with 
momentum k and spin a. The spin-cr fermions have mass 
771ct and dispersion = fc^/2mcr, and V is the volume 
of the system. We take h ~ 1 throughout the paper, and 
consider the mass-balanced case m-f- = = m. All the 
theoretical considerations are for zero temperature (or 
T '^Tp with Tp the Fermi temperature) . The diluteness 
of the system ensures that the range h of the potential 
is much smaller than the typical inter-particle distance 
l/kp, or kpb <C 1, with kp the Fermi momentum of the 
spin-up sea, and therefore the details of the interaction 
potential become irrelevant. Accordingly, without loss 
of generality, one can model the short-ranged interaction 
as a contact interaction, ^(r) = goSir), in combination 
with the standard ultra-violet divergence regularization 
procedure described below. 

The one-body and two-body propagators provide ac- 
cess to information about the 'polaron' and 'molecule' 
channel respectively. The polaron and molecule are two 
distinct objects belonging to different charge sectors. The 



one-body and two-body propagators are discussed in Sec- 
tions HlX] and HTCl The adopted regularization procedure 
for the renormalized interaction is the subject of Section 
HTB] The DiagMC method is introduced in Section HTD] 

A. One-body propagator 

The polaron quasiparticle properties can be extracted 
from the impurity's Green's function defined as 

G,(k,r) = -0(r)(<n5k;W4i(O) l<'> . (3) 

with c^^{t) the annihilation operator in the Heisenberg 
picture, 

The propagator G^(k, r) is written in the momentum 
imaginary-time representation, is a free parameter, N^r 
is the number operator for spin-cr particles, and fj,-f is the 
chemical potential of the spin-up sea. The state 

|<^) = \),\FSiN^)) , (5) 

consists of the spin-down vacuum and the non-interacting 
spin-up Fermi sea. Since we are dealing with an impurity 
spin-down atom, Gj, is only non-zero for times t > 0. The 
ground-state energy and Z-factor can be extracted from 
the Green's function of Eq. Inserting a complete set 
of cigcnstates \^n'^) of the full Hamiltonian ^ for one 
spin-down particle and A^-|- spin-up particles into Eq. 1^ 
yields for k = 

G;(0,r) = -e(r)^|(*^|c^J<^)p 

n 

-^Zp e-(^--'^)^ , (6) 

with Ep the energy of the polaron, En{N^) the energy 
eigenvalues of the Hamiltonian ([2|) and Eps ~ 3 epN^/h 
the energy of the ideal spin- up Fermi gas, with ep = 
kp I (2to) the Fermi energy. 

The difference between the polaronic and molecular 
state is embedded in the factors K^n'^ |cq^|$^^) P in 
Eq. ([H]). For situations where the polaron is a well-defined 
quasiparticle in the ground state l^'g we have Eq. ([T|) 
for the Z-factor and Ep ~ Eo{N^) — Eps- If, on the 
other hand, the ground state 1^*^^) is a dressed molecule 

the overlap (^o'^ Icqj^I^q'^) is zeroii. This is clear from 
the expansion of the molecular state in the number of 
particle- hole excitations, 

i<') = (E'^kcVkt 

^ k 

+ E' ekk'q<_k-k'i4t4't^qt + • • ■ ) ■ (7) 

k,k',q ^ 



The coefficients ^ are variational parameters, and the 
primes indicate that the sums on k, k' and q are re- 
stricted to |k|, |k'| > kp and |q| < kp. Even if a molecule 
is formed in the ground state, the polaron can be a well- 
defined excited state (in the sense of a narrow peak in 
the spectral function), and Zp can be non-zero. 

For vanishing interactions V the impurity Green's 
function of Eq. ^ becomes 



G°(k,r) 



-(ek4.-M)T 



(8) 



The one-body propagator for the spin-up sea is defined 
as 
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FIG. 1: Graphical representation of the Dyson equation and 
the Bethe-Salpeter equation. The free (dressed) one-body 
impurity propagator is denoted by (G^). The E and 11 are 
the one-body and two-body self-energies, respectively. The 
F is the fully dressed interaction, wheres F'^ is the partially 
dressed interaction obtained by summing all the bare ladders 
GjG? (see Eqs. (fT6l) and (fT7|)). 



Gt(k,r) = -{^>^'\Tr[c^^{T)ci,{Q)\ , (9) 

with Tt the time-ordering operator. Without interac- 
tions, one obtains the free propagator 



G°(k,r) = 



-(fk-fF)r0C|J^| 



kp) if r > 
Ikl) if T<0 



Our goal is to calculate the Gj, of Eq. Q to extract 
Ep by means of Eq. This is achieved by summing all 
irreducible one-particle self-energy diagrams with the Di- 
agMC algorithm (which works in momcntum-imaginary- 
time representation). The irreducible self-energy I](k, w) 
in imaginary-frequency representation is obtained after a 
numerical Fourier transform, and inserted into Dyson's 
equation to give G^, 



[G;(k,a.)]-'= [G;(k,^)]-'-S](k,. 



(10) 



with Lo the imaginary frequency. A graphical represen- 
tation of the Dyson equation is shown in the top panel 
of Fig. [TJ As was shown in Refi^, the polaron energy 
Ep and Z-factor Zp can be extracted directly from the 
self-energy I](0,t), 



p + tx> 

Ep ^ dr l](0,r) e^^" 

Jo 



1 -/+°°dT T I](0,t) e(-^=p-^)^ 



(11) 
(12) 



The effective mass of the polaron is evaluated with 
the estimator— 



1/Zp 



1/m + Ba 



with 



Jv2S(k,r) 



k=0 



(13) 



(14) 



which can conveniently be estimated by expanding 
I](k, r) in Legendre polynomials. One obtains 



iv^E(k,r)|..o 
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J^^dkJ:ik,r)(^^-?j ,(1 



5) 



and the integral can be evaluated during the MC simula- 
tion. The upper limit of integration (A) is optimized to 
minimize the statistical noise while avoiding a systematic 
error at too large A. We also used an alternative way by 
calculating the quasiparticle spectrum -E(k) and fitting 
m* via £'(k) = Ep + k'^/ (2m*). 



B. Renormalized interaction 

We introduce the s-wave scattering length a for colli- 
sions between spin-up and spin-down particles. One of 
the advantages of working with Feynman diagrams is that 
one can work directly in the zero-range limit kpb — > (or, 
equivalent ly, A/kp — )■ -f oo with A an ultraviolet momen- 
tum cut-off) while keeping kpa constant. Thereby, the 
ultra-violet physics can be taken into account by means 
of a summation over all Feynman ladder diagrams. 

In momentum-imaginary-frequency representation 
(p, ri), one obtains for the partially dressed interaction 

r"(p, n) = go + go n°(p, n) r°(p, n) , (i6) 

with n*^ the two-particle self-energy consisting of one 
'bare' ladder 

n°(p,i^) = -^E /rf^G?(| + q,..) 



-y 



|q|<A 

2 

0(|p/2 + q|-A:f^) 



X G0(^-q,J7-c.) 



V I -^^^ zfJ — p2/(4m) — jm + fi + ep 



(17) 



where the momentum cutoff A is required to keep 
the sum finite. The bare coupling constant ^(p) = 
/ dve^^^'^V {t) = go can be eliminated in favor of the 
physical scattering length a by using standard scattering 
theory 
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m 
47ra 



k <A 



(18) 



The F°(p, r2) from Eq. ([T5|) can be expressed in terms 
of the s-wave scattering length a, by taking the limit 
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A — ^ +00 and go with a fixed. In this zero-range 
Hmit, one gets 



[rO(p,i]ri = [fO(p,i])]-i-n(p,i7), 



(19) 



with 



J0r ^^^^ - 



r_ 

2m 



(p-q) 

2m 



Here, we have taken the thermodynamic limit (V — > -|-c» 
and N^/V fixed). The integral in Eq. pO| can be evalu- 
ated analytically, and the dressed interaction in vacuum 
is given by 



[f«(p,r!)]-i = - — v/p2-4m(i^]-FA^ + ef), (21) 

47ra OTT 

for ri 7^ or /i < — e_F, and assuming the principal 
branch. For ji < — [ep + l/(ma^)], the Fourier transform 
to imaginary time can be done analytically, producing 



f"(p,r)=- 



47r 
1 



TTT 



ma 



'erfcf — \ 

v V m a 



, (22) 



with erfc(x) the complementary error function. As in 
Ref4, we use r°(p,T) as a partially dressed interaction 
vertex in the diagrammatic series, instead of the bare 
interaction vertex gg. This dressed vertex is calculated 
here in imaginary time representation by performing the 
Fourier transform of Eq. (jl9p numerically. 

In a next step, the interaction vertex will be fully 
dressed by calculating the two-particle self-energy 11 and 
plugging it into the Bethe-Salpeter equation. 



[r(p,f^)] 



[rO(p,i7)]-i-n(p,f]) 



(23) 



A graphical representation of this equation is shown in 
Fig. [TJ The self-energy 11 contains all connected two- 
particle diagrams that are irreducible with respect to cut- 
ting a single F*' propagator. To avoid double counting, 
the diagrams for 11 should not contain any ladders, since 
those have been summed in r° by means of the Eq. 
This rule also holds when summing diagrams for the one- 
body self-energy E, built from free propagators and 



C. Two-body propagator 

Here, we consider the pair annihilation operator, 

A = E^(q)W^qi' (24) 



with (/3(q) the momentum representation of the wave 
function if{r) for relative motion of the two fermions of 
opposite spin. The two-particle propagator is defined as 



G2(k,r) = -%)(<I'o 14(^)^.1(0) , 



(25) 



where we included the fact that the impurity spin- J, atom 
propagates forward in time. Inserting the complete basis 



^n'^'^^) for (A'^^-l- 1) spin-up particles and one spin-down 



(20) particle, gives 



G2(0,r) 



- -^(-)EK*"^^'i^on<ii^ 

n 

X f^-(E„iN-f + l)-EFs{N-f)-,ir-tJ-)r 



(26) 



with Enioi the molecule energy and Z^oi the molecule Z- 
factor. If the molecule is a well-defined quasiparticle in 
the ground state, we have 



(27) 



and E^oi = Eo{N^ + 1) - Efs{N^) - Mr- Note that the 
value of Zmoi depends on the wave function V3(q). The 
functional from of this pair wave function depends on the 
nature of experiment used to probe the molecule. 

In practice, it is easier to calculate the molecule energy 
from the fully dressed interaction F (see Eq. ([^51) ). This 
function is closely related to the pair correlation function, 
namely. 



F(k,r)=.go<5(r)-f.go^(k,r) 50 



(28) 



with 



r{r,r) = -0(r)(<i>^|(M;)(r,r)(*^#|)(O,O)|<1 , 

(29) 

the pair correlation function. The field operators 
^^^(r) = X^k ^ *'^'"^k cr/"^/^ create a spin-cr fermion at 
position r. In Eq. (|29l) the pair of particles is created 
at the same position (which corresponds to <^(q) = 1 
in Eq. ([Ml))- The structure of the fully dressed in- 
teraction F and the two-particle propagator G2 now 
implies that both structures have the same poles (see 
Eq. (pS]) ). Therefore, the exponential tail of the function 
F(k = 0,t) can conveniently be used for estimating the 
molecule energy, rather than the tail of 62(0, r). This is 
equivalent with looking for this pole of the Bethe-Salpeter 
equation ([25)1 . The molecule's energy i?moi is given by the 
parameter fi that satisfies the equation 



[F°(p = o,n = 0)]-^ = H(p = 0, r2 = 0) , 



(30) 



where the left-hand-side is known analytically, and the 
right-hand-side is evaluated with the DiagMC algorithm 
in imaginary-time domain. 



D. Diagrammatic Monte Carlo 

DiagMC evaluates the series of Feynnian diagrams for 
the self-energy in a stochastic way. We deal with both the 
one-body and two-body self-energies. In a first step, the 
self-energy is built from the free propagators G° and the 
partially dressed interaction r° (obtained through sum- 
mation of G°G° ladders; as discussed in Section Hi Bp . 
We will refer to this series as the 'bare series'. Fig. [2] 
shows the one-body and two-body self-energy diagrams 
up to order 3 in the bare scheme. The order of a dia- 
gram is N when there are N dressed interactions F (i.e., 
N boxes) present in the E-diagram, and N — 1 boxes in 
the Il-diagram. Note that the diagrams cannot contain 
ladders since these have been taken into the vertex func- 
tion F". To illustrate the factorial growth with order, the 
number of one-body self-energy diagrams for given order 
A'' < 12 is given in the second column of Table U for the 
bare series. 

In a second step, we will use dressed propagators or 
'bold lines' in the diagrams. Such dressed (skeleton) se- 
ries are evaluated with the Bold DiagMC technique^ii^. 
We consider the case with only dressed Gj, propagators 
while keeping F° of Eqs. ([T^ and (PH]) as renormalized 
interaction, and the case whereby both the one-body 
propagators and interactions are dressed. We will re- 
fer to these skeleton series as 'bold G' and 'bold G-F', 
respectively. In the latter case, the Bold DiagMC algo- 
rithm is constructed as follows: given approximate one- 
body and two-body self-energies E and 11, the Dyson and 
Bethe-Salpeter equation are solved to deliver the one- 
body propagator Gj, and the dressed interaction F (see 
Eqs. pU)) and ([Ml))- In a next step, these are used to 
dress the series for E and 11, which are evaluated stochas- 
tically with DiagMC up to order N^, . This self-consistent 
cycle is repeated until convergence is reached. Fig. [3] 
shows the skeleton (bold G-F) series for the one- and two- 
body self-energies up to order 4. Evidently, when dress- 
ing the lines in the self-energies, one has to keep track 
of two-particle reducibility, and systematically avoid any 
double counting. This typically means that at any order 
N the numbers in the second column of Table U (Bare) 
are an upper limit of the number of diagrams in the third 
and fourth column. At = 2 and = 4, however, the 
number of diagrams increases due to the fact that ladders 
should be included again once Gj, is bold. All the dia- 
grams of Table U are summed explicitly during the (Bold) 
DiagMC simulation. 



III. RESUMMATION AND BOLDIFICATION 

When considering a diagrammatic series, it is natural 
to ask whether there are dominant classes of diagrams. 
Identification of the dominant diagrams potentially al- 
lows one to make good approximations. To address this 
issue, we constructed a histogram counting how many 
times a certain topology is sampled. We consider first 



5 




FIG. 2: Diagrammatic expansion for the one-body self-energy 
E and the two-body self-energy 11. Here, the diagrams are 
built from the bare propagators (thin lines), and the par- 
tially dressed interaction F" (light grey box). All diagrams 
have a 'backbone' structure, since we have a single impurity 
propagating forward in time and interacting with a Fermi sea 
of free particles. 




FIG. 3: Skeleton diagrammatic expansion for the one-body 
and two-body self-energy: the impurity propagator and in- 
teraction lines that appear in the diagram are fully dressed 
solutions of the Dyson equation and the Bethe-Salpeter equa- 
tion (see Fig. 




FIG. 4: The figure contains the two dominant one-body self- 
energy diagrams for N — 6. Imaginary time runs from right 
to left. 
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N 


Bare 


Bold G 


Bold G-r 


1 


1 


1 


1 


2 





1 





3 


2 


2 


1 


4 


6 


7 


2 


5 


34 


34 


13 


6 


210 


206 


74 


7 


1,526 


1,476 


544 


8 


12,558 


12,123 


4,458 


9 


115,618 


111,866 


41,221 


10 


1,177,170 


1,143,554 


421,412 


11 


13,136,102 


12,816,572 


4,722,881 


12 


159,467,022 


156,217,782 


57,553,440 



TABLE I: Factorial increase of the number of Feynman dia- 
grams. At fixed order A'^, the number of one-body self-energy 
diagrams is given for different types of series: the bare series, 
the skeleton series with dressed Gj. (bold G), and the skeleton 
series with dressed G^ and F lines (bold G-F). 



T 




Bare series 
Bold G series 
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FIG. 5; (color online) The polaron energy in units of the 
Fermi energy as a function of the inverse maximum diagram 
order 1/Nt for irreducible self-energy diagrams at unitarity 
l/{kFa) = 0. The red squares show the polaron energy cal- 
culated via Eq. (|lip with self-energy diagrams built from the 
free propagators G'^ and the partially dressed propagator F". 
The blue circles show the results from the bold-G approach. 

the bare series. It turns out that for the one-body self- 
energy, roughly half of the simulation time is spent on 
sampling two diagrams at each order. These two dia- 
grams are shown in Figure 2] for diagram order six. To 
understand why these two diagrams arc dominant at a 
fixed N, we use an argument first made by Hugcnholtzii. 
For the dilute spin-up gas, momentum integration in- 
side the Fermi sea is heavily restricted in phase space 
(momentum integration runs up to the Fermi momen- 
tum kp ~ (N^/^)^^^)- This implies the presence of a 
backward (or hole) spin-up propagator reducing the con- 
tribution of the diagram significantly, while the forward 
(particle) propagator enhances the contribution roughly 



with a factor Jjj,|^j,^dk. As a consequence, diagrams 
with the smallest possible number of hole propagators 
will be dominant. For the self-energy, we see that, at 
fixed order, the minimum number of hole propagators is 
two. Since the number of fermion loops differs by one, 
these two diagrams have opposite sign. Numerically we 
found that the two diagrams almost cancel each other. 
This can be seen in Figure [5l where we show the polaron 
energy Ep as a function of the inverse diagram order cut- 
off N^, for the interaction strength l/{kpa) = 0. For the 
'bare series', we observe a fast convergence due to cancel- 
lation of diagrams. This magic cancellation was referred 
to as 'sign blessing—. At infinite scattering length, such 
near cancellation was also observed by Combescot and 
Giraudi^. They have found that the success of the Chevy 
ansatz at strong coupling can be attributed to a nearly 
perfect destructive interference of the states with more 
than one particle-hole excitation. Combescot and Giraud 
illustrated that an expansion in powers of the hole wave 
vectors converges extremely rapidly at unitarity. In our 
case, the series is organized differently, but at fixed order 
we have exactly the same type of cancellation between 
diagrams with the same number of hole propagators^^. 
Just like in the Combescot- Giraud argument, the cancel- 
lation is exact when the momentum-dependence of the 
hole propagators is neglected. Note that the dominant 
diagrams (see Figure 2]) can also be viewed as three-body 
T-matrix diagrams closed with two hole propagators^. 
This class of diagrams, in which there are at most two 
particle-hole excitations, has been considered previously 
for the polaron proble m^^i^^ . It was shown that they 
exactly reproduce the Skorniakov and Ter-Martirosian 
equationi^ in the BEG limit. In this strong-coupling 
limit, the dominant process is scattering between a dimer 
and a spin-up fermion, which is diagrammatically repre- 
sented by the three-body T-matrix diagrams. Away from 
this limit, the considered class of diagrams turns out to 
give a quantitatively good correction to the lowest order 
result. We find that this is due to their dominance, even 
away from the BEC limit. 

When going towards the BEC side (1/a > 0), the can- 
cellation between dominant diagrams of the type shown 
in Figure m is no longer perfect. Figure [5] shows the po- 
laron energy as a function of for l/{kpa) — 1.333. 
For the bare series, the oscillations prevent one from ex- 
tracting Ep for 1/iV* — >• 0. 

To cure the bad convergence of the bare series for 
l/{kpa) > one can include more diagrams by dress- 
ing the propagators. We start by dressing the spin-down 
propagator lines, while keeping the partially dressed F*^. 
Diagrams reducible with respect to cutting two spin- 
down lines should no longer be sampled, since they are 
included implicitly. For l/{kpa) = the self-energy 
S(k, t) converges in this 'bold G scheme' for N^, > 7. 
Extrapolation to infinite N^, gives the exact S and Gj.. 
Figure [5] includes the polaron energy as function of dia- 
gram order cut-off when the one-body self-energy is built 
with the exact G^. The bare and bold series converge to 
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FIG. 6: (color online) Same plot as Figure O but now consid- 
ering the interaction strength parameter l/ikpa) — 1.333. In 
the bare series, small oscillations unable us to extrapolate to 
infinite diagram order. 

the same energy. Remarkably, the dressed scheme gives 
worse results at low TV*. This indicates that approxima- 
tions based on a few low order diagrams are completely 
uncontrolled, and including more diagrams by dressing 
the lines does not necessarily improve the quality of the 
results. 

For l/{kpa) ~ 1.333, we sec that dressing the impu- 
rity lines helps to get rid of the residual oscillations in 
the bare scheme (see Figure (HI)- One might expect that 
dressing even more, by using a fully dressed F instead 
of r°, might lead to even better convergence. Figure [7] 
shows however that, even for l/(/cFa) = 0, the fully bold 
scries (bold G-F scheme) does not seem to converge {N^, 
is the diagram cut-off for both E and IT, and a Bold Di- 
agMC simulation is done for each iV* ) , in contrast to the 
results of Ref<^. The data for the fully bold simulation 
of Refj^ was obtained by using the exact G and F (i.e. 
extrapolated to the — >■ oo limit with resummation fac- 
tors). They were not obtained with a self-consistent sim- 
ulation, which explains the difference. Moreover, data is 
not shown above N^, = 7, where oscillations do occur. In 
order to understand why the series no longer converges, 
we introduce an intermediate scheme (which we call bold 
G-F^): the self-energy is built from the fully converged 
Gj. and a partially dressed interaction F^, built from sum- 
ming the ladders GiG^. The result is also shown in Fig- 
ure [71 and we again observe convergence to the same 
answer as in Figure [51 The key difference between both 
schemes is that in the bold G-F^ scheme, both dominant 
diagrams shown in Fig. Estill explicitly contribute to the 
self-energy, whereas in the fully bold scheme the upper 
dominant diagram becomes reducible and is taken into 
account self-consistently. This means that the balance 
of cancellation between diagrams is broken, and a sin- 
gle dominant diagram keeps contributing at each order. 
So, it turns out that dressing the diagrams as much as 
possible is not always a good idea. In this respect, our 



FIG. 7: (color online) The polaron energy in units of the 
Fermi energy as a function of the inverse maximum diagram 
order 1/Nt, for irreducible self-energy diagrams at unitarity 
l/ikpo) = 0. The blue circles show results from a fully Bold 
DiagMC simulation: the self-energy diagrams are built from 
fully dressed one-particle propagators G4. and two-particle 
propagators F up to self-energy diagram order A'^,. The red 
squares show the polaron energy calculated with diagrams 
built with the exact and a partially dressed interaction 
containing the sum of all G^G^ ladders. 

findings disagree with Refi^. 

A second method to cure the bad convergence of the 
bare series on the BEG side, is to employ series resum- 
mation techniques. We will use the Abelian resummation 
techniques'^ which have been used when calculating the 
equation of state of the unitary gas with Bold DiagMG^S. 
This resummation technique works as follows. One starts 
from a series f{x) = ^„ that has a finite radius of 
convergence R > 0. The idea is to sum the series at some 
point xq outside of the radius R by analytically contin- 
uing the function /. This provides a good procedure for 
summing the divergent series in the sense that it respects 
basic operations (sum, multiplication and derivative) and 
that it preserves distinctness^^. It is well-known that 
with analytic continuation, one can encounter problems 
with the existence and/or uniqueness of the solutior*2i. 
However, one can formally define a domain called the 
'Mittag-Leffier star' where the function can be analyti- 
cally continued along straight lines [0, xq]. Note that this 
star will always contain the disk of convergence. It can 
be showrti^ that for each point xq of the Mittag-Leffier 
star, the limit 

lim y dnx^e-'^" , (31) 

n 

with Xn — n log(n) for n > and Aq = 0, exists and is 
equal to the analytic continuation of / to the point xq. 
Note that within the disk of convergence the procedure 
works equally well, and can improve the rate of conver- 
gence. We apply the Abelian resummation technique to 
the expansion of the self-energy S and 11. As the analytic 
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FIG. 8: (color online) Illustration of the Abelian resumma- 
tion technique for the geometric series. We evaluate /^{xo) = 
Xge"'^" for xo = —3. and various choices of A„. The 
value of the analytically continued function 1/(1 — sq) is re- 
trieved for e — >■ 0"'' . 



FIG. 9: (color online) Abelian resummation of the bare series 
of one-body self-energy diagrams at l/ikpo) = 1.333. The 
polaron energy Ep/ep is extracted in the limit e = 0^ for 
different choices of A„. 



structure of S and 11 is unknown, it is currently impossi- 
ble to determine whether there is a finite radius of con- 
vergence and whether we are in the Mittag-LefHer star. 
In practice, we apply different resummation techniques 
(i.e., different functions A„), and test the uniqueness of 
the result. 

We use following A„: (i) Lindelof 1: A„ = n log(n) for 
n > and Aq = 0; (ii) Lindelof 2: A„ = (n— 1) log(n— 1) 
for n > 1 and Aq = Ai = 0; (iii) Gauss 1: A„ = ri^ 
for n > 0; (iv) Gauss 2: A„ = (n — 1)^ for n > 1 and 
Aq = 0; (iv) Gauss 3: A,i = (n — 2)^ for n > 2 and 
Aq = Ai = 0. Before applying these resummation tech- 
niques to our diagrammatic series, we illustrate its power 
with an example for the geometric series. Figure |8] shows 
the sums /e(xo = —3) = XqC"'^" for a few choices of 
A„. At small e. the computation of /e is no longer fea- 
sible due to finite computer precision. By extrapolating 
to e = 0, we indeed find 1/(1 — xq) with high precision. 
The Lindelof curve gives a slighly less accurate extrapo- 
lation because it suppresses high order contributions in a 
much smoother fashion than the Gaussian resummation. 
When applying these techniques to our diagrammatic se- 
ries, it is the growth of the statistical error bars (due 
to factorial complexity) that prevents us from going to 
very small values of e. Figure [9] shows the polaron energy 
calculated with the resummed self-energy as a function 
of the control parameter e for l/(fci?a) = 1.333. The 
polaron energy Ep can be extracted with high accuracy. 
The major source of error bar stems from the uncertainty 
in the extrapolation. 

Histogramming the different topologies of the two- 
body self-energies H revealed a dominant diagram at 
each order. This diagram is shown in Figure 1101 Again 
it shows a three-body T-matrix structure that is closed 
with a single spin-up hole propagator. Upon increasing 
the diagram order up to 20, we observe a steady growth 



FIG. 10: At fixed order, there is one dominant diagram for 
the two-body self-energy. Here, we draw this diagram at order 
six. 



in the contribution of this diagram. This is illustrated 
in Figure llll where we plot the n-th order contribution 
n„ to the two-body self-energy as a function of imag- 
inary time r for external momentum zero. Figure [T^] 
illustrates that we can nonetheless get accurate values 
for the molecule energy -Emoi by using different Abelian 
resummation techniques and extrapolating to e = 0+. 
Again, the Gaussian resummation methods allow one to 
reach very small values of e. The quoted error bars are 
rather conservative as we include the extrapolated results 
obtained with all choices for A„. 

We also tested the resummability of the fully bold se- 
ries (bold G - bold F scheme), since this was used in the 
Bold DiagMC method for determining the equation of 
state of the unitary gas^. When applying the Abelian 
resummation techniques and extrapolate to e = 0"*" at 
l/(/ci?a) = 0, the correct polaron energy is retrieved. 
This constitutes an independent check for the resumma- 
tion of the skeleton series. 



IV. QUASIPARTICLE PROPERTIES 

As an independent cross-check of Refi^, which uses al- 
ternate ways of resumming the diagrammatic series, we 
calculate the ground-state energies of the polaron and 
molecule. Figure [T3] shows these energies shifted by the 
vacuum molecule energy E}, = — l/(ma^) in units of the 
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FIG. 11: (color online) The two-particle self-energy 11 at ex- 
ternal momentum zero as a function of imaginary time for 
kpa = 1. The n„ are the contribution of all n-th order di- 
agrams, and here shown for various n. As n increases, n„ 
keeps on growing. We work in units kp = ^, m = 1, h = 1 
and yL/ep = —3.2. The noise in the curves indicates the mag- 
nitude of the statistical error. 



FIG. 13: (color online) The extracted polaron and molecule 
energy as a function of the interaction strength l/(fcFa). The 
energies are expressed as {E — Ef,)/eF, with Ei, = — l/(ma^) 
the molecule energy in vacuum. The FN-DMC results are 
from Refi^, the variational results for the polaron from 
Ref.— and for the molecule from Ref.— . DiagMC results by 
Prokof'ev and Svistuno\^ are also shown. 




0.2 0.4 0.6 0.8 1 

e 



FIG. 12: (color online) Abelian resummation of the bare series 
of two-particle self-energy diagrams aX kpa = 1. The molecule 
energy -Emoi/tF is extracted in the limit e = 0^ for different 
choices of A„. 



Fermi energy tp- A selection of the polaron and molecule 
energies is also given in Table |lll We find the transi- 
tion point at (kpa)c = 1.15(3), in agreement with Ref4. 
Close to the transition point, we find polaron energies 
that differ about 1% with the polaron energies of Rcfi^, 
which, we believe, is due to a small systematic error in the 
lowest order diagram in Ref.—. The variational energies 
obtained from a wave-function ansatz for the polaron^ 
and the moleculeii are very close to the Monte Carlo 
results. Note that Chevy's variational ansatz for the 
polaron state is completely equivalent with the non-self- 
consistent T- matrix approximatior^i^ which is exactly our 
bare series at A^* = 1. Fixed node-difussion Monte Carlo 



(FN-DMC) results are also in good agreement with the 
DiagMC data. For l/{kFa) = 2 it seems that system- 
atic errors on the FN-DMC results were underestimated, 
since FN-DMC should in principle give an upper bound 
to the true ground-state energy. 

Figure [T3] shows the effective mass of the polaron 
as calculated with DiagMC. We compare with the 
ENS experiment^ at unitarity, DiagMC calculations by 
Prokof 'ev-SvistunoA*^, FN-DMC— li, a variational calcula- 
tion up to two particle-hole excitationsi^, and the first or- 
der (A^* = 1) result in the bare scheme and the fully bold 
G-T scheme. The experimental effective mass, which is 
in perfect agreement with DiagMC^, was extracted from 
the low frequency breathing modes, and in particular 
the Fermi polaron breathing mode. The lowest order 
bare calculation, also known as T-matrix approximation, 
is equivalent to the Chevy ansatz, while the lowest or- 
der bold calculation corresponds to the self-consistent T- 
matrix approximation. These results show that including 
only single particle-hole pair excitions does not lead to 
accurate results for the effective mass, while the varia- 
tional calculation based on diagrams taking into account 
at most two particle-hole pairs excitations agrees with 
the DiagMC resultsi^. 

Experimental and theoretical quasiparticle residues are 
shown in Fig. 1151 To create and probe polarons, the MIT 
experiment^ starts from a cloud of ®Li atoms with most 
atoms occupying the lowest hyperfinc state |1) (spin- up), 
and about 2% of the atoms occyping the hyperfinc state 
|3) (spin-down) in the degenerate regime T « 0.14 Tp 
with Tp the Fermi temperature. A broad Feshbach res- 
onance is used to enhance the scattering between atoms 
in states |1) and |3). Radio- frequency (rf) spectra of the 
spin-up and spin-down components arc measured. The 
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FIG. 14: (color online) The effective mass m* of the po- 
laron in units of the bare mass m as a function of the in- 
teraction parameter l/{kFa). Our DiagMC results (open cir- 
cles) are shown together with DiagMC results by Prokof'ev 
and Svistuno\«^ (filled blue circles), FN-DMC results^ii (black 
triangles), ENS experimenli^ (red square) and a variational 
calculation up to two particle-hole excitationsi^ (solid black 
line). We also show m, calculated from the lowest order 
self-energy diagram (i.e., Nt = 1) for the bare series (solid 
grey line) and for the fully bold G-F series (solid blue line), 
which are equivalent to the non-self-consistent and the self- 
consistent T-matrix approximation, respectively. 

atoms arc transferred to a third empty state with very 
weak final-state interactions. Therefore, the measured 
transition rate / can be connected with the impurity's 
spectral function in linear response theory^i^ 

/(wl) oc ^ 7iF(ek " - wl) Pi(k, ek - /i - ujl) , (32) 

k 

with lol the frequency of the rf photons and np^x) = 
1/(1-1- e^^) the Fermi distribution. Note that the spec- 
tral fmiction depends on the temperature. Density in- 
homogeneities are taken care of through tomographic 
reconstruction^. At sufficiently weak attractions, the 
Fermi polaron is observed as a narrow peak in the im- 
purity spectrum that is not matched by the broad envi- 
ronment spectrum. The peak position gives the polaron 
energy Ep, and was found to be in perfect agreement 
with the DiagMC results of Ref<^. The polaron Z-factor 
was measured by determining the ratio of the area under 
the impurity peak that is not matched by the environ- 
ment, and the total area under the impurity's spectrum. 
The experimental Z-factor from Refj^ is shown in Figure 
[TSl together with the Z-factor calculated from Chevy's 
ansataiii^, the fully self-consistent result in lowest order 
(iV* = 1) and our DiagMC simulation. DiagMC data for 
the factor is also given in Table HIl 

The results obtained via DiagMC simulation agree ex- 
tremely well with Chevy's variational ansatz. This is 
very surprising in the strongly interacting regime where 
Zp is significantly smaller than one. Here, one would ex- 



pect multiple particle-hole excitations to be important 
since the overlap with the non-interacting wave-function 
is small. Remarkably, including just single particle-hole 
excitations on top of the Fermi sea produces almost the 
exact Zp . When the lowest order diagram is calculated in 
a fully self-consistent way, however, the agreement with 
DiagMC is less good. This hints at the fact that the 
almost perfect agreement with Chevy's ansatz (i.e., the 
lowest order bare result) is rather accidental. 

The Z-factors computed with Chevy's ansatz and Di- 
agMC both exceed the measured ones. It was sug- 
gested in Refjii that the disagreement between the ex- 
periment and the Chevy ansatz is an artefact of Chevy's 
expansion being restricted to one particle-one hole ex- 
citations. As the DiagMC technique includes multiple 
particle-multiple hole excitations and agrees very well 
with Chevy's ansatz, we see that this is not the case. 
However, since the measured Zp might only give a lower 
bound^, theory and experiment might not be in disagree- 
ment. 

The measured polaron Z-factor vanishes beyond a 
critical interaction strength. Ignoring issues related to 
metastability, once the two-body bound state becomes 
energetically favorable, all polarons disappear and the 
polaron peak vanishes. In the experiment T/Tp = 
0.14(3), and finite-temperature effects are thus expected 
to become important. Indeed, close to {kpa)c = 1.15(3) 
the energy difference between the molecule and polaron 
state is of the order of 0.1 Tp (Fig. fT5|) . Therefore, one ex- 
pects that T = calculations underestimate the critical 
\/{kpa) measured at T « 0.1 Tp. Similarly, the mea- 
sured l/{kpa)c can be interpreted as an upper bound 
for the r = situation. On the other hand, due to de- 
pletion of the experimental spectrum, the measured Z 
might only give a lower bound, which means that the 
experimentally determined critical l/(kpa) might be un- 
derestimated. These uncertainties might explain why the 
critical 1/ {kpo) in the experiment is lower than the value 
obtained with DiagMC for a single impurity. Fixed-node 
Monte-Carlo simulations for a finite density of impurities, 
on the other hand, predict phase separation before the 
systems even reaches the polaron-to-molecule transition^, 
and the vanishing Z-factor might be a manifestation of 
this phase separation. 



V. CONCLUSIONS 

We have considered the Fermi-polaron system in three 
dimensions, in which a single spin-down impurity is 
strongly coupled to a non-interacting Fermi sea of spin- 
up particles. Although this system contains strongly in- 
teracting fermions, it can be solved with the Diagram- 
matic Monte-Carlo method. This method is based on the 
stochastic evaluation of a series of Feynman diagrams. 
To extract ground-state properties, one has to overcome 
a factorial complexity due to the increase of the number 
of diagrams. Nonetheless, extrapolation to infinite dia- 
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FIG. 15: (color online) The polaron quasiparticle residue Zp 
as a function of the interaction parameter l/(fc_Fa). DiagMC 
results (open circles) are compared with variational ansataii 
(black solid line), the fully bold G-T series at A^* = 1 or self- 
consistent T-matrix approximation (solid blue line) and the 
MIT experimenli^ (blue circles). 



l/ikpa) 


Ep/Ep 


Emol/Ep 


Zp 


-1.8 


-0.1793(1) 




0.9727(4) 


-1.6 


-0.1961(1) 




0.9665(5) 


-1.4 


-0.2159(2) 




0.9590(3) 


-1.2 


-0.2393(2) 




0.9502(3) 


-1.0 


-0.2687(2) 




0.9376(4) 


-0.8 


-0.3052(2) 




0.9209(5) 


-0.6 


-0.3526(2) 




0.8978(8) 


-0.4 


-0.4141(2) 




0.8670(10) 


-0.2 


-0.4976(2) 




0.8237(15) 


0.0 


-0.615(1) 




0.7586(27) 


0.2 


-0.782(1) 




0.6720(42) 


0.4 


-1.028(2) 




0.5672(28) 


0.6 


-1.385(2) 


-1.180(13) 


0.4410(32) 


0.8 


-1.880(2) 


-1.830(8) 


0.3258(58) 


1.0 


-2.540(3) 


-2.618(6) 


0.2283(70) 


1.2 


-3.372(4) 


-3.554(6) 


0.1559(69) 


1.4 


-4.373(5) 


-4.633(5) 


0.1102(68) 


1.6 


-5.554(8) 


-5.867(6) 


0.0771(58) 


1.8 


-6.889(12) 


-7.251(5) 


0.0578(35) 



TABLE II: Selection of DiagMC data for the polaron energy 
Ep, molecule energy -Bmoi and polaron residue Zp for several 
values of the interaction strength parameter l/{kFa). 



gram order becomes possible when the diagrams cancel 
each other better than the factorial increase in number. 
At interaction strength l/ikpa) = 0, we find such per- 
fect cancellation (within our statistical errors). When 
considering the series built on bare propagators on the 
BEC side, however, oscillations with diagram order re- 
main and prevent a controlled extrapolation to the in- 
finite diagram order. We have followed two strategies 
around this problem: the first is to consider skeleton se- 
ries (built on dressed propagators), and the second the 
use of resumniation techniques. Though dressed series 
can be evaluated to higher orders, we have found that 
in some cases dressing can destroy a favorable cancella- 
tion of diagrams. For all interaction strengths we found 
that the (skeleton) series of the one-body and two-body 
self-energy is rcsummable by means of Abelian resum- 
mation. Bare series, skeleton series and resummcd series 
give robust answers in their respective region of applica- 
bility (i.e. where the infinite diagram order extrapolation 
is controlled). 

We have identified classes of dominant diagrams for the 
one-body and two-body self-energy in the crossover re- 
gion of strong interaction. The dominant diagrams turn 
out to be the leading processes of the strong-coupling 
limit: scattering between a dimer and a spin-up fermion, 
which is diagrammatically represented by the three-body 
T-matrix diagrams. Including just these dominant dia- 
grams gives a quantitatively good correction to the lowest 
order result, even away from the strong-coupling limit. 

We have shown that not only the polaron and molecule 
energies agree very well with a variational ansatz from 
weak to strong attraction, but also the polaron residue or 
Z-factor. Though this agreement must be due to strong 
cancellation of diagrams, we only observed convergence 
for the bare series at l/{kFa) = 0. A full explanation 
for the success of the variational ansatz is still missing, 
and it is therefore unclear in which cases the ansatz is 
appropriate. 
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